/* 

CRIME IS IN THE AIR: THE CONTEMPORANEOUS RELATIONSHIP BETWEEN AIR POLLUTION AND CRIME

Malvina Bondy, Sefi Roth, Lutz Sager

October 2019
	
This program contains all Stata routines used to generate the Figure and Tables in the manuscript as published.	

Contact: Sefi Roth (S.J.Roth@lse.ac.uk)

*/		



*--------------------------------------------------------------------------
* Preamble
*--------------------------------------------------------------------------
	* Loading the Data
		clear all
		set more off
		set maxvar 8000
		set matsize 10000
		* working directory
			cd ".../Crime and Air Quality/"
		* set temp directory for intermediate results
			global TEMP ".../temp"	
		* load data
			use "03_CrimeIsInTheAir.dta"

	* Defining the Set of Control Variables & Instruments
		* Control Variables: Rainfall, Temperature (bins), Relative Humidity, Relathive Humidity squared, Interaction of Temperature and Relative Humidity, Squared Interaction of Relative Humidity and Temperature, Windspeed, Police Deployment, Unemployment, Tube Activity, Dummy for period after July 2005 Terror Attack (post-treatment in Draca et al. 2011), Dummy for Wards in Focus of Police Redeployment following the Attack (treated in Draca et al. 2011), Interaction of the previous two (treated wards during treated period in Draca et al. 2011)
			global controls "LS3dw_rainfall i.LS3dw_temperature_bin LS3dw_relhumidity LS3dw_rh2 LS3dw_temp_rh LS3dw_temp2_rh2 LS3dw_windspeed lpolicepop un tube i.afterattack i.policy6 i.treat1_policy6"
		* Wind Instruments: 	3 out of 4 wind direction quadrants, separately for each of 5 regions
			global instwind "LS3dw_WIND_DIRDUM_1_901 LS3dw_WIND_DIRDUM_91_1801 LS3dw_WIND_DIRDUM_271_3601 LS3dw_WIND_DIRDUM_1_902 LS3dw_WIND_DIRDUM_91_1802 LS3dw_WIND_DIRDUM_271_3602 LS3dw_WIND_DIRDUM_1_903 LS3dw_WIND_DIRDUM_91_1803 LS3dw_WIND_DIRDUM_271_3603 LS3dw_WIND_DIRDUM_1_904 LS3dw_WIND_DIRDUM_91_1804 LS3dw_WIND_DIRDUM_271_3604 LS3dw_WIND_DIRDUM_1_905 LS3dw_WIND_DIRDUM_91_1805 LS3dw_WIND_DIRDUM_271_3605"
		* Inversion Instrument:	The continuous temperature difference between station-level temperature and low-level atmospheric temperature
			global instinv "NASA_D_DIFF_station"
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------


*--------------------------------------------------------------------------
* Table 1: Descriptive Statistics
*--------------------------------------------------------------------------

* 3 iterations: (1) full sample, (2) wards with below median house prices, (3) wards with above median house price
	preserve
		quietly reg crime LS3dw_AQI_10, robust
		keep if e(sample)
		quietly outreg2 using "Table_1.xls", replace sum(log) eqkeep(N mean sd) label keep(WARDavg_HP_2004_2005 lcrime crime_perpop pop_ward LS3dw_AQI LS3dw_AQIuk LS3dw_PM10 LS3dw_CO LS3dw_NO2 LS3dw_SO2 LS3dw_O3 LS3dw_temperature LS3dw_rainfall LS3dw_relhumidity LS3dw_windspeed lpolicepop un tube)
	restore
	preserve
		quietly reg crime LS3dw_AQI_10, robust
		keep if e(sample)
		keep if WARD_q4_HP <=1
		quietly outreg2 using "Table_1.xls", append sum(log) eqkeep(N mean sd) label keep(WARDavg_HP_2004_2005 lcrime crime_perpop pop_ward LS3dw_AQI LS3dw_AQIuk LS3dw_PM10 LS3dw_CO LS3dw_NO2 LS3dw_SO2 LS3dw_O3 LS3dw_temperature LS3dw_rainfall LS3dw_relhumidity LS3dw_windspeed lpolicepop un tube)
	restore
		preserve
		quietly reg crime LS3dw_AQI_10, robust
		keep if e(sample)
		keep if WARD_q4_HP >=2
		quietly outreg2 using "Table_1.xls", append sum(log) eqkeep(N mean sd) label keep(WARDavg_HP_2004_2005 lcrime crime_perpop pop_ward LS3dw_AQI LS3dw_AQIuk LS3dw_PM10 LS3dw_CO LS3dw_NO2 LS3dw_SO2 LS3dw_O3 LS3dw_temperature LS3dw_rainfall LS3dw_relhumidity LS3dw_windspeed lpolicepop un tube)
	restore
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
	
*--------------------------------------------------------------------------
* Table 2: Pooled and Fixed Effect Models of Air Pollution’s Impact on Crime
*--------------------------------------------------------------------------

	* (1) PPML: No controls (no fixed effects, simply absorbing a constant one=1)
		ppmlhdfe crime LS3dw_AQI_10, absorb(one) cluster(ward_code yearmonth) offset(log_ONSpop)

			* generate normalised coefficient and elasticity estimate
				quietly mat b = e(b)
				quietly scalar coef = b[1,1]
				quietly scalar pseudor2 = e(r2_p)
				quietly sum LS3dw_AQI_10 if e(sample)
				quietly scalar stdev = r(sd)
				quietly scalar avg = r(mean)
				quietly scalar coefnorm = exp(coef * stdev)-1
				quietly scalar elast = coef / (1/avg)
			* write to file
				quietly outreg2 using "Table_2.xls", replace bdec(3) sdec(4) ctitle (POIS) keep(LS3dw_AQI_10) addstat(Pseudo R2, pseudor2, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, NO, Ward FE, NO, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Year-Months)

	* (2) PPML: Controls
		ppmlhdfe crime LS3dw_AQI_10 $controls, absorb(one) cluster(ward_code yearmonth) offset(log_ONSpop)

			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly scalar pseudor2 = e(r2_p)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_2.xls", append bdec(3) sdec(4) ctitle (POIS) keep(LS3dw_AQI_10 $controls) addstat(Pseudo R2, pseudor2, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, NO, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Year-Months)
		
	* (3) PPML: Controls + Ward Fixed Effects
		ppmlhdfe crime LS3dw_AQI_10 $controls, absorb(ward_code) cluster(ward_code yearmonth) offset(log_ONSpop)

			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly scalar pseudor2 = e(r2_p)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_2.xls", append bdec(3) sdec(4) ctitle (POIS) keep(LS3dw_AQI_10 $controls) addstat(Pseudo R2, pseudor2, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Year-Months)

	* (4) PPML: Controls + Ward FE + Day-of-Week FE
		ppmlhdfe crime LS3dw_AQI_10 $controls, absorb(ward_code dow) cluster(ward_code yearmonth) offset(log_ONSpop)

			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly scalar pseudor2 = e(r2_p)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_2.xls", append bdec(3) sdec(4) ctitle (POIS) keep(LS3dw_AQI_10 $controls) addstat(Pseudo R2, pseudor2, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, NO, SE cluster, Wards/Year-Months)

	* (5) PPML: Controls + Ward FE + DOW FE + Year-by-Month FE
		ppmlhdfe crime LS3dw_AQI_10 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)

			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly scalar pseudor2 = e(r2_p)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_2.xls", append bdec(3) sdec(4) ctitle (POIS) keep(LS3dw_AQI_10 $controls) addstat(Pseudo R2, pseudor2, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)

	* (6) OLS (crime per capita as dependent variable):  Controls + Ward FE + DOW FE + Year-by-Month FE
		quietly reghdfe crime_perpop LS3dw_AQI_10 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth)

			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly sum crime_perpop if e(sample)
			quietly scalar crimeavg = r(mean)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = (coef/crimeavg) * stdev
			quietly scalar elast = (coef/crimeavg) / (1/avg)
			quietly outreg2 using "Table_2.xls", append bdec(3) sdec(4) ctitle (OLS (crime_perpop)) keep(LS3dw_AQI_10 $controls) addstat(Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
	
*--------------------------------------------------------------------------
* Table 3: Air Pollution's Impact on Crime – Non-linear Models
*--------------------------------------------------------------------------

	preserve
			* split AQI levels into 5 groups with similar bracket size
			capture drop LS3dw_AQI_g5
			quietly egen LS3dw_AQI_g5=cut(LS3dw_AQI), at(0,20,25,30,35,105)
			label define AQIgroup 0 "Group 1 (<20)" 20 "Group 2 (20 - 25)" 25 "Group 3 (25 - 30)" 30 "Group 4 (30 - 35)" 35 "Group 5 (>35)"
			label values LS3dw_AQI_g5 AQIgroup
			*(1)
			quietly ppmlhdfe crime ib0.LS3dw_AQI_g5, absorb(one) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly outreg2 using "Table_3.xls", replace bdec(3) sdec(4) ctitle (OLS) keep(20.LS3dw_AQI_g5 25.LS3dw_AQI_g5 30.LS3dw_AQI_g5 35.LS3dw_AQI_g5) addtext(Ward FE, NO, Weather controls, NO, Other controls, NO, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Month-Year)
			*(2)
			quietly ppmlhdfe crime ib0.LS3dw_AQI_g5 $controls, absorb(one) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly outreg2 using "Table_3.xls", append bdec(3) sdec(4) ctitle (OLS) keep(20.LS3dw_AQI_g5 25.LS3dw_AQI_g5 30.LS3dw_AQI_g5 35.LS3dw_AQI_g5) addtext(Ward FE, NO, Weather controls, YES, Other controls, YES, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Month-Year)
			*(3)		
			quietly ppmlhdfe crime ib0.LS3dw_AQI_g5, absorb(ward_code) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly outreg2 using "Table_3.xls", append bdec(3) sdec(4) ctitle (OLS) keep(20.LS3dw_AQI_g5 25.LS3dw_AQI_g5 30.LS3dw_AQI_g5 35.LS3dw_AQI_g5) addtext(Ward FE, YES, Weather controls, NO, Other controls, NO, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Month-Year)
			*(4)			
			quietly ppmlhdfe crime ib0.LS3dw_AQI_g5 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly outreg2 using "Table_3.xls", append bdec(3) sdec(4) ctitle (OLS) keep(20.LS3dw_AQI_g5 25.LS3dw_AQI_g5 30.LS3dw_AQI_g5 35.LS3dw_AQI_g5) addtext(Ward FE, YES, Weather controls, YES, Other controls, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Month-Year)
	restore		
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
			
*--------------------------------------------------------------------------
* Table 4: Air Pollution's Impact on Crime – Instrumental Variable Models using Wind Direction
*--------------------------------------------------------------------------

	* (1) PPML (not instrumented): Controls + Ward FE + DOW FE + Year-by-Month FE
		ppmlhdfe crime LS3dw_AQI_10 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly scalar pseudor2 = e(r2_p)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_4.xls", replace bdec(3) sdec(4) ctitle (PPML) keep(LS3dw_AQI_10 $controls) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)
			
	* (2) PPML (not instrumented): Controls + Ward FE + DOW FE + Year-by-Month FE []sub-sample that has instruments]	
			quietly capture drop AQI_10_inst
			quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow i.yearmonth, absorb(ward_code) cluster(ward_code yearmonth)
			quietly predict AQI_10_inst
		ppmlhdfe crime LS3dw_AQI_10 $controls if AQI_10_inst !=., absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly outreg2 using "Table_4.xls", append bdec(3) sdec(4) ctitle (PPML) keep(LS3dw_AQI_10 $controls) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)

	* (3) PPML (instrumented): Controls

		* FIRST STAGE:
		quietly capture drop AQI_10_inst
		quietly reghdfe LS3dw_AQI_10 $instwind $controls, absorb(one) cluster(ward_code yearmonth)
			quietly test $instwind
			quietly scalar ftest = r(F)
			quietly outreg2 using "Table_A2.xls", append bdec(3) sdec(4) ctitle (OLS) keep($instwind) addstat(Instrument F-test, ftest) addtext(controls, YES, Ward FE, NO, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Year-Months)
			quietly predict AQI_10_inst

		* SECOND STAGE:
		ppmlhdfe crime AQI_10_inst $controls, absorb(one) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_4.xls", append bdec(3) sdec(4) ctitle (2SLS PPML) keep(AQI_10_inst $controls) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, NO, DOW FE, NO, Year-Month FE, NO, SE cluster, Wards/Year-Months)

	* (4) PPML (instrumented): Controls + Ward FE + DOW FE

		* FIRST STAGE:
		quietly capture drop AQI_10_inst
		quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow, absorb(ward_code) cluster(ward_code yearmonth)
			quietly test $instwind
			quietly scalar ftest = r(F)
			quietly outreg2 using "Table_A2.xls", append bdec(3) sdec(4) ctitle (OLS) keep($instwind) addstat(Instrument F-test, ftest) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, NO, SE cluster, Wards/Year-Months)
			quietly predict AQI_10_inst

		* SECOND STAGE:
		ppmlhdfe crime AQI_10_inst $controls, absorb(ward_code dow) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_4.xls", append bdec(3) sdec(4) ctitle (2SLS PPML) keep(AQI_10_inst $controls) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, NO, SE cluster, Wards/Year-Months)

	* (5) PPML (instrumented): Controls + Ward FE + DOW FE + Year-by-Month FE

		* FIRST STAGE:
		quietly capture drop AQI_10_inst
		quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow i.yearmonth, absorb(ward_code) cluster(ward_code yearmonth)
			quietly test $instwind
			quietly scalar ftest = r(F)
			quietly outreg2 using "Table_A2.xls", append bdec(3) sdec(4) ctitle (OLS) keep($instwind) addstat(Instrument F-test, ftest) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)
			quietly predict AQI_10_inst

		* SECOND STAGE:
		ppmlhdfe crime AQI_10_inst $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = exp(coef * stdev)-1
			quietly scalar elast = coef / (1/avg)
			quietly outreg2 using "Table_4.xls", append bdec(3) sdec(4) ctitle (2SLS PPML) keep(AQI_10_inst $controls) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)
		
	* (6) 2SLS (instrumented per capita crime):  Controls + Ward FE + DOW FE + Year-by-Month FE

		* FIRST STAGE:
		quietly capture drop AQI_10_inst
		quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow i.yearmonth, absorb(ward_code) cluster(ward_code yearmonth)
			quietly test $instwind
			quietly scalar ftest = r(F)
			quietly outreg2 using "Table_A2.xls", append bdec(3) sdec(4) ctitle (OLS) keep($instwind) addstat(Instrument F-test, ftest) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)
			quietly predict AQI_10_inst

		* SECOND STAGE:
		reghdfe crime_perpop AQI_10_inst $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth)	
			quietly mat b = e(b)
			quietly scalar coef = b[1,1]
			quietly sum crime_perpop if e(sample)
			quietly scalar crimeavg = r(mean)
			quietly sum LS3dw_AQI_10 if e(sample)
			quietly scalar stdev = r(sd)
			quietly scalar avg = r(mean)
			quietly scalar coefnorm = (coef/crimeavg) * stdev
			quietly scalar elast = (coef/crimeavg) / (1/avg)
			quietly outreg2 using "Table_4.xls", append bdec(3) sdec(4) ctitle (2SLS (crime_perpop)) keep(AQI_10_inst LS3dw_rainfall i.LS3dw_temperature_bin LS3dw_relhumidity LS3dw_rh2 LS3dw_temp_rh LS3dw_temp2_rh2 LS3dw_windspeed lpolicepop un tube i.afterattack i.policy6  i.treat1_policy6) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Ward FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Wards/Year-Months)
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
		
*--------------------------------------------------------------------------
* Table 5: Air Pollution's Impact on Crime – Air Pollution’s Impact on Crime – Different Levels of Aggregation
*--------------------------------------------------------------------------
						 
	* PANEL A: BOROUGH LEVEL 	
		preserve 
					bysort ward_code : egen ward_totalcrimes = total(crime)
					collapse (rawsum) ONSpersons crime* police  (mean) policepop LS3dw_AQI_10 NASA_D_INV NASA_temp_D2 NASA_temp_D1 LS3dw_rainfall LS3dw_temperature LS3dw_relhumidity LS3dw_windspeed un tube afterattack policy6 treat1 treat1_policy6 dow month yearmonth LS3dw_WIND* /* [aweight = ward_totalcrimes]*/, by(borough LS_date)

					capture drop lcrime_no0
					gen lcrime_no0 = log(crime)
					gen lpolicepop = log(policepop)
					quietly egen LS3dw_temperature_bin=cut(LS3dw_temperature), group(5)
					quietly gen LS3dw_rh2 = LS3dw_relhumidity ^ 2
					quietly gen LS3dw_temp_rh = LS3dw_temperature * LS3dw_relhumidity
					quietly gen LS3dw_temp2_rh2 = LS3dw_temp_rh ^ 2
					quietly gen LS3dw_temp2 = LS3dw_temperature ^ 2
					gen NASA_temp_D2_C = NASA_temp_D2 - 273.15
					capture drop NASA_D_DIFF_station
					gen NASA_D_DIFF_station = NASA_temp_D2_C - LS3dw_temperature
					capture gen log_ONSpop = log(ONSpersons)
					capture gen one = 1
					drop if NASA_D_DIFF_station ==.

			* (A1) PPML (not instrumented): Controls + Borough FE + DOW FE + Year-by-Month FE
					ppmlhdfe crime LS3dw_AQI_10 $controls, absorb(borough dow yearmonth) cluster(yearmonth) offset(log_ONSpop)
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", replace bdec(3) sdec(4) ctitle (PPML) keep(LS3dw_AQI_10 ) addstat(Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)
						 
			* (A2) PPML (instrumented using wind): Controls + Borough FE + DOW FE + Year-by-Month FE
					global Binstwind "LS3dw_WIND_DIRDUM_1_901 LS3dw_WIND_DIRDUM_91_1801 LS3dw_WIND_DIRDUM_181_2701 LS3dw_WIND_DIRDUM_1_902 LS3dw_WIND_DIRDUM_91_1802 LS3dw_WIND_DIRDUM_181_2702 LS3dw_WIND_DIRDUM_1_903 LS3dw_WIND_DIRDUM_91_1803 LS3dw_WIND_DIRDUM_181_2703 LS3dw_WIND_DIRDUM_1_904 LS3dw_WIND_DIRDUM_91_1804 LS3dw_WIND_DIRDUM_181_2704 LS3dw_WIND_DIRDUM_1_905 LS3dw_WIND_DIRDUM_91_1805 LS3dw_WIND_DIRDUM_181_2705"

					* FIRST STAGE:
					quietly capture drop AQI_10_inst
					quietly reghdfe LS3dw_AQI_10 $Binstwind $controls i.dow i.yearmonth, absorb(borough) cluster(yearmonth)
						quietly test $Binstwind
						quietly scalar ftest = r(F)
					quietly predict AQI_10_inst
					* SECOND STAGE:
					ppmlhdfe crime AQI_10_inst $controls, absorb(borough dow yearmonth) cluster(yearmonth) offset(log_ONSpop)
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (2SLS PPMLwind}wind) keep(AQI_10_inst ) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)

			* (A3) PPML (instrumented using thermal inversions): Controls + Borough FE + DOW FE + Year-by-Month FE
					global Binstinv "NASA_D_DIFF_station"

					* FIRST STAGE:
					quietly capture drop AQI_10_inst
					quietly reghdfe LS3dw_AQI_10 $Binstinv $controls i.dow i.yearmonth, absorb(borough) cluster(yearmonth)
						quietly test $Binstinv
						quietly scalar ftest = r(F)
					quietly predict AQI_10_inst
					* SECOND STAGE:
					ppmlhdfe crime AQI_10_inst $controls, absorb(borough dow yearmonth) cluster(yearmonth) offset(log_ONSpop)
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (2SLS PPMLinv) keep(AQI_10_inst ) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)

			* (A4) PPML (instrumented using both wind and inversions): Controls + Borough FE + DOW FE + Year-by-Month FE
					global Binstboth "NASA_D_DIFF_station LS3dw_WIND_DIRDUM_1_901 LS3dw_WIND_DIRDUM_91_1801 LS3dw_WIND_DIRDUM_181_2701 LS3dw_WIND_DIRDUM_1_902 LS3dw_WIND_DIRDUM_91_1802 LS3dw_WIND_DIRDUM_181_2702 LS3dw_WIND_DIRDUM_1_903 LS3dw_WIND_DIRDUM_91_1803 LS3dw_WIND_DIRDUM_181_2703 LS3dw_WIND_DIRDUM_1_904 LS3dw_WIND_DIRDUM_91_1804 LS3dw_WIND_DIRDUM_181_2704 LS3dw_WIND_DIRDUM_1_905 LS3dw_WIND_DIRDUM_91_1805 LS3dw_WIND_DIRDUM_181_2705"

					* FIRST STAGE:
					quietly capture drop AQI_10_inst
					quietly reghdfe LS3dw_AQI_10 $Binstboth $controls i.dow i.yearmonth, absorb(borough) cluster(yearmonth)
						quietly test $Binstboth
						quietly scalar ftest = r(F)
					quietly predict AQI_10_inst
					* SECOND STAGE:
					ppmlhdfe crime AQI_10_inst $controls, absorb(borough dow yearmonth) cluster(yearmonth) offset(log_ONSpop) d
						quietly capture drop temp_*
						quietly predict temp_Yhat
						quietly gen temp_uhat = crime - temp_Yhat
						quietly reg temp_uhat $Binstboth
						quietly scalar hansenj = e(r2)*e(N)
						quietly scalar hansenp = 1-chi2(3, hansenj)				
					ppmlhdfe crime AQI_10_inst $controls, absorb(borough dow yearmonth) cluster(yearmonth) offset(log_ONSpop)
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (2SLS PPMLboth) keep(AQI_10_inst ) addstat(Overid J, hansenj, Overid p, hansenp, K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, YES, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)
		restore
							
	* PANEL B: LONDON LEVEL 	
		preserve 
					bysort ward_code : egen ward_totalcrimes = total(crime)
					collapse (rawsum) ONSpersons crime* police  (mean) policepop LS3dw_AQI_10 NASA_D_INV NASA_temp_D2 NASA_temp_D1 LS3dw_rainfall LS3dw_temperature LS3dw_relhumidity LS3dw_windspeed un tube afterattack policy6 treat1 treat1_policy6 dow month yearmonth LS3dw_WIND* /* [aweight = ward_totalcrimes]*/, by(LS_date)

					capture drop lcrime_no0
					gen lcrime_no0 = log(crime)
					gen lpolicepop = log(policepop)
					quietly egen LS3dw_temperature_bin=cut(LS3dw_temperature), group(5)
					quietly gen LS3dw_rh2 = LS3dw_relhumidity ^ 2
					quietly gen LS3dw_temp_rh = LS3dw_temperature * LS3dw_relhumidity
					quietly gen LS3dw_temp2_rh2 = LS3dw_temp_rh ^ 2
					quietly gen LS3dw_temp2 = LS3dw_temperature ^ 2
					gen NASA_temp_D2_C = NASA_temp_D2 - 273.15
					capture drop NASA_D_DIFF_station
					gen NASA_D_DIFF_station = NASA_temp_D2_C - LS3dw_temperature
					capture gen log_ONSpop = log(ONSpersons)

					drop if NASA_D_DIFF_station ==.

					global loncontrols "LS3dw_rainfall i.LS3dw_temperature_bin LS3dw_relhumidity LS3dw_rh2 LS3dw_temp_rh LS3dw_temp2_rh2 LS3dw_windspeed lpolicepop un tube i.afterattack i.policy6"
					
			* (B1) PPML (not instrumented): Controls + DOW FE + Year-by-Month FE
					ppmlhdfe crime LS3dw_AQI_10 $loncontrols, absorb( dow yearmonth) cluster( yearmonth) offset(log_ONSpop)					
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (PPML) keep(LS3dw_AQI_10 ) addstat(Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, NO, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)
						 
			* (B2) PPML (instrumented using wind): Controls + DOW FE + Year-by-Month FE
					global Binstwind "LS3dw_WIND_DIR_1_90 LS3dw_WIND_DIR_91_180 LS3dw_WIND_DIR_271_360"

				* FIRST STAGE:
					quietly capture drop AQI_10_inst
					quietly reghdfe LS3dw_AQI_10 $Binstwind $loncontrols , absorb(dow yearmonth) cluster(yearmonth) 
						quietly test $Binstwind
						quietly scalar ftest = r(F)
					quietly predict AQI_10_inst

				* SECOND STAGE:
					ppmlhdfe crime AQI_10_inst $loncontrols , absorb(dow yearmonth) cluster(yearmonth) offset(log_ONSpop) 					
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (2SLS POISwind) keep(AQI_10_inst ) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, NO, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)

			* (B3) PPML (instrumented using thermal inversions): Controls + DOW FE + Year-by-Month FE
					global Binstinv "NASA_D_DIFF_station"

				* FIRST STAGE:
					quietly capture drop AQI_10_inst
					quietly reghdfe LS3dw_AQI_10 $Binstinv $loncontrols , absorb(dow yearmonth) cluster(yearmonth) 
						quietly test $Binstinv
						quietly scalar ftest = r(F)
					quietly predict AQI_10_inst

				* SECOND STAGE:
					ppmlhdfe crime AQI_10_inst $loncontrols, absorb(dow yearmonth) cluster(yearmonth) offset(log_ONSpop) 
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (2SLS POISinv) keep(AQI_10_inst ) addstat(K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, NO, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)

			* (B4) PPML (instrumented using both wind and inversions): Controls + Borough FE + DOW FE + Year-by-Month FE
					global Binstboth "NASA_D_DIFF_station LS3dw_WIND_DIR_1_90 LS3dw_WIND_DIR_91_180 LS3dw_WIND_DIR_271_360"

				* FIRST STAGE:
					quietly capture drop AQI_10_inst
					quietly reghdfe LS3dw_AQI_10 $Binstboth $loncontrols , absorb(dow yearmonth) cluster(yearmonth) 
						quietly test $Binstboth
						quietly scalar ftest = r(F)
					quietly predict AQI_10_inst
					
				* SECOND STAGE:
					ppmlhdfe crime AQI_10_inst $loncontrols, absorb(dow yearmonth) cluster(yearmonth) offset(log_ONSpop) d	
						quietly capture drop temp_*
						quietly predict temp_Yhat
						quietly gen temp_uhat = crime - temp_Yhat
						quietly reg temp_uhat $Binstboth
						quietly scalar hansenj = e(r2)*e(N)
						quietly scalar hansenp = 1-chi2(3, hansenj)				
						* alternatively: scalar hansenj = 4*e(F)
					ppmlhdfe crime AQI_10_inst $loncontrols, absorb(dow yearmonth) cluster(yearmonth) offset(log_ONSpop)				
						quietly mat b = e(b)
						quietly scalar coef = b[1,1]
						quietly sum LS3dw_AQI_10 if e(sample)
						quietly scalar stdev = r(sd)
						quietly scalar avg = r(mean)
						quietly scalar coefnorm = exp(coef * stdev)-1
						quietly scalar elast = coef / (1/avg)
						quietly outreg2 using "Table_5.xls", append bdec(3) sdec(4) ctitle (2SLS POISboth) keep(AQI_10_inst ) addstat(Overid J, hansenj, Overid p, hansenp, K-P F-test, ftest, Beta normalised, coefnorm, Elasticity, elast) addtext(controls, YES, Borough FE, NO, DOW FE, YES, Year-Month FE, YES, SE cluster, Year-Months)
			restore
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
				 
*--------------------------------------------------------------------------
* Figure 1: Geographic coverage – London wards, AURN/MIDAS stations
*--------------------------------------------------------------------------
	* Map drawn using QGIS software
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
		
*--------------------------------------------------------------------------
* Figure 2: AQI and Crime – Density plots
*--------------------------------------------------------------------------
	quietly hist LS3dw_AQI, scheme(s1color) lcolor(black) color(navy) bin(100)
	graph export "Figure_2a.png", replace
	quietly hist crime if crime <=50, scheme(s1color) lcolor(black) color(navy) bin(50)
	graph export "Figure_2b.png", replace
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------

*--------------------------------------------------------------------------
* Figure 3: Residual AQI and Crime - Binned scatterplot
*--------------------------------------------------------------------------
	preserve
		set more off

		global weather "LS3dw_rainfall i.LS3dw_temperature_bin LS3dw_relhumidity LS3dw_rh2 LS3dw_temp_rh LS3dw_temp2_rh2 LS3dw_windspeed "

		xtset ward_code LS_date
		xtreg lcrime_no0 $weather ,fe cluster(ward_code)
		predict predcrime
		gen residcrime=lcrime_no0-predcrime
		xtreg LS3dw_AQI $weather ,fe cluster(ward_code)
		predict predaqi
		gen residaqi = LS3dw_AQI - predaqi
		reg residcrime residaqi
		gen tresidaqi = (_b[residaqi]/_se[residaqi])
		gen bresidaqi=_b[residaqi]
		gen presidaqi = (2 * ttail(e(df_r), abs(_b[residaqi]/_se[residaqi])))
		global bresidaqi=round(bresidaqi,.001)

		 global tresidaqi=round(tresidaqi,.01)
		 global presidaqi=round(presidaqi,.0001)
		 global mychar="="
		 if presidaqi<.001{
		 global presidaqi="<.001"
		 global mychar=""
		 }
		 
		drop if residcrime==.
		drop if residaqi==.
		gen n=_N
		global n=n
		gen myrand=uniform()
		egen aqicat=cut(residaqi),at(-40(3)100)
		drop if aqicat<-35
		drop if aqicat>50

		collapse (mean) residcrime,by(aqicat) fast
		lowess residcrime aqicat ,gen(residsmooth)
		graph twoway (scatter residcrime aqicat ,msymbol(oh) msize(*1.5) mcolor(black)) (line residsmooth aqicat, lwidth(*1.5) lcolor(black) xtitle("{bf:Residual AQI}",height(6)) ytitle("{bf:Residual Log of Crime}",height(3)) yscale(titlegap(*+10)) legend(off)), scheme(s1color)
		graph export "Figure_3.png", replace
	restore
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------

*--------------------------------------------------------------------------
* Figure 4: House prices - London wards
*--------------------------------------------------------------------------
	* Map drawn using QGIS software
	* Long-run averages of ward-level house prices as follows
	preserve
		collapse (mean)  WARDavg_HP_2004_2005, by(gsscode)
		rename gsscode GSS_CODE
		export delimited using "Figure_4_Data.csv", replace
	restore
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------

*--------------------------------------------------------------------------
* Figure 5: The Effect of Air Pollution by House Price Quantile
*--------------------------------------------------------------------------
	* Separate regression for each HP quintile. Looks similar with one regression with interaction terms.
	preserve

		quietly egen groups = cut(WARDavg_HP_2004_2005), group(5)
		replace groups = groups + 1
		matrix median = J(5, 3, .)
		matrix rown median = 1 2 3 4 5 
		matrix mean = J(5, 3, .)
		matrix rown mean = 1 2 3 4 5 
		forvalues i=1(1)5 {
			quietly centile LS3dw_AQI if groups == `i'
			matrix median[`i',1] = r(c_1), r(lb_1), r(ub_1)
			quietly mean LS3dw_AQI if groups == `i'
			matrix mtemp`i' = e(b)
			matrix mean[`i',1] = mtemp`i'[1,1]
			matrix mtump`i' = e(V)
			matrix mean[`i',2] = mtemp`i'[1,1] - 1.96* mtump`i'[1,1]^(1/2)
			matrix mean[`i',3] = mtemp`i'[1,1] + 1.96* mtump`i'[1,1]^(1/2)
		}

		matrix beta = J(5,3,.)
		matrix rown beta = 1 2 3 4 5 

		forvalues i=1(1)5 {

			* FIRST STAGE:
			quietly capture drop AQI_10_inst
			quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow i.yearmonth, absorb(ward_code) cluster(ward_code yearmonth)
			quietly predict AQI_10_inst	
		
			* SECOND STAGE:
			quietly ppmlhdfe crime AQI_10_inst $controls  if groups == `i', absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
				matrix temp`i' = e(b)
				local coef = temp`i'[1,1]
				matrix beta[`i',1] = `coef'
				matrix tump`i' = e(V)
				matrix beta[`i',2] = temp`i'[1,1] - 1.96* tump`i'[1,1]^(1/2)
				matrix beta[`i',3] = temp`i'[1,1] + 1.96* tump`i'[1,1]^(1/2)
		}
		* Generating the graph from stored results
		coefplot (matrix(beta[,1]), ci((beta[,2] beta[,3])) label(Coefficient estimate (by quantile))) (matrix(mean[,1]), msymbol(D) ci((mean[,2] mean[,3])) axis(2) label(Mean AQI (by quantile))), scheme(s1color) vertical xtitle("{bf:House price decile}") ytitle("{bf:Coefficient estimate per 10 AQI}") ytitle("{bf: Average AQI 2004-05}", axis(2)) yline(0)
		graph export "Figure_5.png", replace

		
	restore
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------
	
*--------------------------------------------------------------------------
* Figure 6: The Effect of Air Pollution by Major Crime Type
*--------------------------------------------------------------------------
		* FIRST STAGE:
		quietly capture drop AQI_10_inst
		quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow i.yearmonth, absorb(ward_code) cluster(ward_code yearmonth)
		quietly predict AQI_10_inst	
		* generating AQI variables with names to indicate model for graph
			capture drop yy*
			capture gen yy1=AQI_10_inst
			capture gen yy2=AQI_10_inst
			capture gen yy3=AQI_10_inst
			capture gen yy4=AQI_10_inst
			capture gen yy5=AQI_10_inst
			capture gen yy6=AQI_10_inst

		* SECOND STAGES:	
		quietly ppmlhdfe crime yy1 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
		estimates store IV2All_Crime
		quietly ppmlhdfe theft yy2 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
		estimates store IV2Theft
		quietly ppmlhdfe violence yy3 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
		estimates store IV2Violence_incsex
		quietly ppmlhdfe damage yy4 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
		estimates store IV2Damage
		quietly ppmlhdfe burglary yy5 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
		estimates store IV2Burglary
		quietly ppmlhdfe robbery yy6 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop)
		estimates store IV2Robbery
		
		coefplot (IV2All_Crime, label(off) offset(0)) (IV2Theft, label(off) offset(0)) (IV2Violence_incsex, label(off) offset(0)) (IV2Damage, label(off) offset(0)) (IV2Burglary, label(off) offset(0)) (IV2Robbery, label(off) offset(0)), keep(yy*) nolabels xtitle() coeflabels(yy1="All Crimes" yy2="Theft" yy3="Violence (incl. sexual offences)" yy4="Criminal Damage" yy5="Burglary" yy6="Robbery", notick labsize(small) wrap(12) labcolor(purple) labgap(1)) omitted scheme(s1color) vertical ytitle("{bf:Coefficient estimate per 10 AQI}") label(off) legend(off) yline(0) 
		graph export "Figure_6.png", replace
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------

*--------------------------------------------------------------------------
* Figure 7: The Effect of Air Pollution by Crime Sub-Type
*--------------------------------------------------------------------------
	* FIRST STAGE:
		quietly capture drop AQI_10_inst
		quietly reghdfe LS3dw_AQI_10 $instwind $controls i.dow i.yearmonth, absorb(ward_code) cluster(ward_code yearmonth)
		quietly predict AQI_10_inst

	* SECOND STAGES:	
		* generating AQI variables with names to indicate model for graph
		capture drop inst1 inst2 inst3 inst4 inst5 inst6 inst7 inst8 inst9 inst10 inst11
		gen inst1 = AQI_10_inst
		gen inst2 = AQI_10_inst
		gen inst3 = AQI_10_inst
		gen inst4 = AQI_10_inst
		gen inst5 = AQI_10_inst
		gen inst6 = AQI_10_inst
		gen inst7 = AQI_10_inst
		gen inst8 = AQI_10_inst
		gen inst9 = AQI_10_inst
		gen inst10 = AQI_10_inst
		gen inst11 = AQI_10_inst 		

		quietly ppmlhdfe burglary_residential inst1 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store Burglary_Residential
		quietly ppmlhdfe burglary_commercial inst2 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store Burglary_Commercial
		quietly ppmlhdfe theft_snatchshoppickpocketcycle inst3 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store theft_other
		quietly ppmlhdfe theft_MV inst4 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store theft_MV
		quietly ppmlhdfe violent_incsex_serious inst5 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store Violent_incsex_Serious
		quietly ppmlhdfe violent_incsex_minor inst6 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store Violent_incsex_Minor
		quietly ppmlhdfe damage_mv inst7 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store damage_mv
		quietly ppmlhdfe damage_other inst8 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store damage_other
		quietly ppmlhdfe crime_4 dd5 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store robbery1
		quietly ppmlhdfe crime_26 inst9 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store robbery2
		quietly ppmlhdfe crime_8 inst10 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store damage_building1
		quietly ppmlhdfe crime_9 inst11 $controls, absorb(ward_code dow yearmonth) cluster(ward_code yearmonth) offset(log_ONSpop) level(90)
		estimates store damage_building2

	coefplot (Violent_incsex_Serious, label(off) offset(0)) (Violent_incsex_Minor, label(off) offset(0)) (damage_building2, label(off) offset(0)) (damage_building1, label(off) offset(0)) (damage_mv, label(off) offset(0)) (theft_other, label(off) offset(0)) (theft_MV, label(off) offset(0)) (Burglary_Residential, label(off) offset(0)) (Burglary_Commercial, label(off) offset(0)) (robbery2, label(off) offset(0))    , keep(inst1 inst2 inst3 inst4 inst5 inst6 inst7 inst8 inst9 inst10 inst11) groups(inst5 inst6="{bf:Violent crimes (inc. sexual)}" inst1 inst2="{bf:Burglary}" inst3 inst4="{bf:Theft}"  inst9="{bf:Robbery}" inst7 inst8 inst11 inst10 ="{bf:Criminal Damage}", labsize(vsmall)) nolabels xtitle() coeflabels(inst5="Serious harm (ABH, GBH, murder, rape)" inst6="Other (Common assault, harassment, other sexual)" inst1="Residential" inst2="Nonresid."  inst3="Snatches, shoplift, pickpocket, cycles, other" inst4="Theft of motor vehicle" inst11="Residential" inst10="Nonresid." inst7="Motor vehicles" inst8="Other" inst9="Personal Property", notick labsize(tiny) wrap(12) labcolor(purple) labgap(2)) omitted scheme(s1color) vertical ytitle("{bf:Coefficient estimate per 10 AQI}") label(off) legend(off) yline(0) xline(8.5) 
	graph export "Figure_7.png", replace
*--------------------------------------------------------------------------
* END
*--------------------------------------------------------------------------


